#include<iostream>
#include<stdlib.h>
#include<iomanip>
#include<ctime>
#include<windows.h>
#include<pthread.h>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX2
using namespace std;
const int n = 2000;
float A[n][n];
float** B;
float** C;
float** D;
void init(float A[n][n])
{
    srand((unsigned)time(NULL));
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
        {
            A[i][j] = 0;
        }
        A[i][i] = 1.0;
        for (int j = i + 1; j < n; j++)
        {
            A[i][j] = rand() % 10000;
        }
    }
    for (int k = 0; k < n; k++)
    {
        for (int i = k + 1; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A[i][j] += A[k][j];
                A[i][j] = float(((int)A[i][j]) % 10000);
            }
        }
    }
}
void display(float** A)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cout << setw(11) << A[i][j] << " ";
        }
        cout << endl;
    }
}
void Gauss_normal_clock(float** A)
{
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}
void clone(float A[n][n], float** B)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
}
typedef struct
{
    int k;
    int t_id;
}threadParam_t;
void* pthread_funct_normal(void* param)
{
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;
    int t_id = p->t_id;
    int i = k + t_id + 1;
    for (int j = k + 1; j < n; j++)
    {
        C[i][j] = C[i][j] - C[i][k] * C[k][j];
    }
    C[i][k] = 0;
    pthread_exit(NULL);
}
void* pthread_funct_sse(void* param)
{
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;
    int t_id = p->t_id;
    int i = k + t_id + 1;
    __m128 vaik = _mm_set1_ps(B[i][k]);
    int j;
    for (j = k + 1; j + 4 <= n; j += 4)
    {
        __m128 vakj = _mm_loadu_ps(&B[k][j]);
        __m128 vaij = _mm_loadu_ps(&B[i][j]);
        __m128 vx = _mm_mul_ps(vakj, vaik);
        vaij = _mm_sub_ps(vaij, vx);
        _mm_storeu_ps(&B[i][j], vaij);
    }
    for (; j < n; j++)
    {
        B[i][j] = B[i][j] - B[k][j] * B[i][k];
    }
    B[i][k] = 0;
    pthread_exit(NULL);
}
void Gauss_pthread_dynamic_normal()
{
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
        {
            C[k][j] = C[k][j] / C[k][k];
        }
        C[k][k] = 1.0;
        int worker_count = n - k - 1;
        pthread_t* handles = (pthread_t*)malloc(worker_count * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(worker_count * sizeof(threadParam_t));
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_create(&handles[t_id], NULL, pthread_funct_normal, (void*)&param[t_id]);
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_join(handles[t_id], NULL);
        }
    }
}
void Gauss_pthread_dynamic_sse()
{
    for (int k = 0; k < n; k++)
    {
        __m128 vt = _mm_set1_ps(B[k][k]);
        int j;
        for (j = k + 1; j + 4 <= n; j += 4)
        {
            __m128 va = _mm_loadu_ps(&B[k][j]);
            va = _mm_div_ps(va, vt);
            _mm_storeu_ps(&B[k][j], va);
        }
        for (; j < n; j++)
        {
            B[k][j] = B[k][j] / B[k][k];
        }
        B[k][k] = 1.0;

        int worker_count = n - k - 1;
        pthread_t* handles = (pthread_t*)malloc(worker_count * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(worker_count * sizeof(threadParam_t));
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_create(&handles[t_id], NULL, pthread_funct_sse, (void*)&param[t_id]);
        }
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            pthread_join(handles[t_id], NULL);
        }
    }
}
int main()
{
    init(A);
    B = new float* [n];
    C = new float* [n];
    D = new float* [n];
    for (int i = 0; i < n; i++)
    {
        B[i] = new float[n];
    }
    for (int i = 0; i < n; i++)
    {
        C[i] = new float[n];
    }
    for (int i = 0; i < n; i++)
    {
        D[i] = new float[n];
    }

    double time1_sum = 0;
    int count1 = 0;
    for (int l = 0; l <10 ; l++)
    {
        double time1 = 0.0;
        clone(A, D);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_normal_clock(D);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time1 = (tail - head) * 1000 / freq;
        if (time1 != 0)count1++;
        time1_sum += time1;
    }
    cout << "普通算法总运行时间：" << time1_sum / count1 << "ms" << endl;

    double time2_sum = 0;
    int count2 = 0;
    for (int l = 0; l < 10; l++)
    {
        double time2 = 0.0;
        clone(A, C);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_pthread_dynamic_normal();
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time2 = (tail - head) * 1000 / freq;
        if (time2 != 0)count2++;
        time2_sum += time2;
    }
    cout << "普通算法pthread动态线程优化总运行时间：" << time2_sum / count2 << "ms" << endl;

    double time3_sum = 0;
    int count3 = 0;
    for (int l = 0; l < 10; l++)
    {
        double time3 = 0.0;
        clone(A, B);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_pthread_dynamic_sse();
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time3 = (tail - head) * 1000 / freq;
        if (time3 != 0)count3++;
        time3_sum += time3;
    }
    cout << "普通算法SSE+pthread优化总运行时间：" << time3_sum / count3 << "ms" << endl;
}